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We reconsider the radial Saffman- Taylor instability, when a fluid injected from a point 
source displaces another fluid with a higher viscosity in a Hele-Shaw cell, where 
the fluids are confined between two neighboring flat plates. The advancing fluid 
front is unstable and forms fingers along the circumference. The so-called Brinkman 
equations is used to describe the flow field, which also takes into account viscous 
stresses in the plane and not only viscous stresses due to the confining plates like the 
Darcy equation. The dispersion relation agrees better with the experimental results 
than the classical linear stability analysis of radial fingering in Hele-Shaw cells that 
uses Darcy's law as a model for the fluid motion. 
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I. INTRODUCTION 



Viscous fingering, also called Saffman- Taylor- instability, is considered as an archetype 
of pattern forming instability (see Couder 20002 for an insightful review). It has also been 
widely studied in the context of industrial research, such as petroleum extraction in partic- 
ular. The phenomenon belongs to the broad family of instabilities of growth in Laplacian 
fields, which includes solidification, aggregation, etc... It was first studied by Saffman and 
Taylor 1 in 1958, who observed the formation of patterns upon injection of a fluid into a 
channel filled with another more viscous fluid. Saffman and Taylor- studied the formation 
of fingers in shallow rectangular channels and observed that the formation of fingers was de- 
pendent on the ratio of channel height to channel width. The extension to radial geometry, 
depicted in figure [Jji, dates back to Bataille 3 , Wilson^ and Paterson^. 




FIG. 1. Sketch of an Hele-Shaw cell with an inviscid liquid, injected amid a viscous liquid. A 
close-up shows schematically the isocontours of the pressure field. 

Viscous and capillary forces govern the mechanism of this instability. The viscosity 
difference drives the formation of fingers, because emerging fingers enhance the pressure 
gradients acting in the tip region differentially in each fluid, as illustrated by iso-contours 
of the pressure field in figure [lb, where the inner fluid is assumed inviscid and at constant 

2 



pressure. Since the interface moves with a velocity proportional to the pressure gradient, 
the feedback loop is closed. In presence of a positive radial gradient of viscosity across the 
interface, any minute initial displacement will turn unstable, as it induces a pressure gradient 
difference that further accelerates a protruding finger, which in turn becomes steeper and 
steeper as it continues to develop. 

This instability is damped by surface tension, which acts to minimize the interface area 
and opposes the formation of curved fingers. One key parameter that determines the number 
of fingers is therefore the ratio between viscous effects, represented by the product of the 
interface velocity U and the viscosity fi, against capillary effects represented by the surface 
tension 7. This ratio is called the capillary number Ca = — ^. 

A second influential parameter is the aspect ratio k = yVlR/h (the origin of the factor 
y/l2 is explained section II). The viscous drag is mainly created by the velocity gradients in 
the small direction h, whereas the capillary forces act on fingers formed along the circum- 
ference whose length scales with the radius R (see figure [T^)- 

Considering a flat front separating fluid 1 to the left from fluid 2 to the right and propa- 
gating along the x coordinate at a velocity V in a rectangular cell of span b (in the direction 
y) and width h (in the direction z), the above physical reasoning has been made rigorous 
by deriving a dispersion relation relating the growth-rate a of the instability to the wave 
number of the disturbance k, assuming a normal mode expansion exp (iax + at). Chuoke 
et al. 7 could demonstrate that the growth-rate was simply given by: 



This formula perfectly reflects the physical picture above and shows the existence of a most 
unstable wavelength, resulting from the competition of a destabilization of the viscosity 
difference and stabilization by surface tension. A similar situation holds in the circular 
injection shown by Wilson^, who determined the wave number of the fastest growing 

perturbation: 



where, for brevity, the injected fluid is considered inviscid and unperturbed. Maxworthy^ 
showed experimentally that this prediction works well for low capillary numbers but does 
not reproduce the behavior at elevated capillary numbers (Ca > 1CT 2 ). 
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Both Chuoke et aF&Z- and Wilson's^ predictions are based on the Darcy equation, which 
is the classical Hele-Shaw limit of the Stokes equation for large aspect ratios. In this limit 
only viscous terms in the thin direction are retained. The flow in the thin direction is 
approximated by a Poiseuille profile and the equations can be depth-averaged in the shallow 
direction. This leads to the Darcy equation, which is a two-dimensional potential flow, 
where the pressure represents the potential. One paradox of this model was pointed out by 
Dai and Shelley-, who argued that for zero surface tension the fingers would form infinitely 
sharp cusps. 

Paterson^ assumed that for zero surface tension flows, for instance with two miscible 
fluids, the interface is influenced by the full three-dimensional stress tensor. He used a 
potential flow for which he derived the 3D stress tensor at the fluid-fluid interface and 
looked for the perturbation wavelength that minimizes the dissipated energy, thereby using 
the property that low Reynolds number flows minimize the dissipation. He found the finger 
size A to be solely dependent on the cell height h in the form A ~ 4.0h. This matches 
quite well the relation Homsy^ deduced from Maxworthy's results: A ~ 5. Oh. A similar 
dissipation minimization approach was recently followed by Boudaoud et a/.—, both for 2D 
Darcy and 2D Stokes suction flows. 

Kim et al.— developed a dispersion relation using viscous potential flow where the normal 
stress boundary condition includes not only pressure but also viscous normal stresses. Their 
result improved the agreement with Maxworthy's experiment even at capillary numbers of 
order one. In this work the Brinkman equation is used as a model for confined flows. The 
Brinkman equation retains in-plane stress terms and is therefore of second order, which al- 
lows boundary conditions for the normal and tangential velocities at the fluid-fluid interface. 
Following the same steps as Paterson^ in his article for radial fingering using Darcy equa- 
tions, we are able to show that the finger prediction in the range from low to high capillary 
numbers is improved. 

In a recent work Logvinov et al.— used the Brinkman equation to describe the Saffman- 
Taylor instability of miscible fluids, hence with zero surface tension, in rectangular channels 
and demonstrate the influence of viscous effects in the flow plane. They found a wavelength 
dependence that is about A ~ 2.5h, which is still by a factor 2 too small compared to the 
experiment. Earlier, Fernandez et alM as well as Carles et al— have also used the Brinkman 
equation to analyze Rayleigh Taylor instability in Hele-Shaw cells. 
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The next section [TT] presents first the Brinkman model and thereafter the stability analysis 
in section III II In section IIVI the results are shown and compare to existing theoretical and 
experimental results. The paper ends by section |VJ where conclusions and perspectives are 
drawn. 



II. MATHEMATICAL FORMULATION OF THE PROBLEM 

The flow geometry under consideration is a Hele-Shaw cell, a space confined between 
two plates, whose distance is generally in the order of a millimeter. The cell is filled with 
a viscous liquid and in the center is a second fluid of smaller viscosity with a uniformly 
circular interface. As the injection of more fluid in the center is started the front advances 
and the instability is triggered. This configuration is schematically depicted in figure QJi. 

The Brinkman equation is used as a model for the fluid motion in this study. Starting 
from the 3D Stokes equation, it is assumed that the velocity component in the thin direction 
is zero and the velocity profiles in the plane direction are Poiseuille profiles (see Boos and 
Thess 1996— and Gallaire et al. 2013^ for more details). Under this condition one can 
integrate in the small direction and obtains depth-averaged equations that combine viscous 
effects in the plane like in the 2D Stokes equation and the viscous friction due to confinement 
from top and bottom wall like in the Darcy equation. The parameter related to the aspect 
ratio that comes from depth averaging is given as k = s/VlR/h, where R is the initial radius 
and h the height of the cell. This gives y/l2 because the derivative of the parabolic velocity 
profile in the thin direction is gives 8R 2 /h 2 , since averaged quantities are used one multiplies 
with the maximum to mean velocity, which is 1.5 for a parabola. 

In constant flux Q radial injection, the interface evolves as R = \J~QtpK. For the non- 
dimensionalisation, using the undisturbed instantaneous radius of the interface p, the flow 
rate Q, the viscosity of the outer fluid This defines the non-dimensional viscosity r\i in 
eq.Q as r\\ = pi/p2 = V an d V2 = 1- This further defines following characteristic velocity 
and pressure: 

U=®- and P Q ^ 



2nR 2ttR 2 ' 

The non-dimensional Brinkman equation then reads 



rjj(Auj — k 2 Uj) — Vpj — with V • = 0. (3) 
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Herein the velocity components of vector u., and pressure pj are only two-dimensional, as 
all the operators. The bottom index corresponds to j — 1 for the inner fluid and j = 2 for 
the outer flow. Given the geometry of the problem it is convenient to write the Brinkman 
equation in polar coordinates: 

(d 2 Ui ldu-j ld 2 Uj iij , o 2d 2 v,j\ dp* 
^{-^- + -^7 + -^--^- k \ ~ ^) ~^ = ° ( 4 ) 

/d\ | 1 dVj | 1 d\ Vj ^ | 2d 2 Uj \ 1 dpj = Q 

V dr 2 r dr r 2 d6 2 r 2 r 2 d6 2 ) r d6 



where Uj and v j denote respectively the radial and azimuthal velocity components. With 

■ dipj _ dipj 
rd9 ' dr 



help of the stream function, Uj = — ^r), a more compact form is obtained by taking 



the curl of eq. (jl]) and eq. ([5]) : 

fld^d d 2 \ / Id r d d 2 fc 2V If _ / 6 x 
\rdr dr r 2 d9 2 / \r9r dr r 2 d6 2 J 



The kinematic boundary conditions at the interface are equality of normal and tangential 
velocity, and normal velocity equal to interface velocity: 

Wi ■ n = V^ 2 • n = p , V^i • t = W 2 • t. (7) 

Since the surface stress component in the shallow direction is averaged out, only the 
radial stress a rr and the in-plane tangential stress a r $ are matched. The dynamic boundary 
conditions are continuous tangential stress: 

t • a\ ■ n = t - <72 • n, (8) 

And jump in normal stresses according to Laplace's law: 

n.a 1 .n-n.a 2 .n = — + (9) 

The corrected boundary conditions of Park and Homsy 17 for an advancing menisc of a 
non wetting fluid is applied, where the in-plane curvature is corrected by 7r/4 

III. LINEAR PERTURBATION THEORY 

The stability of the unperturbed interface is analyzed. The unperturbed interface position 
is p(to) = 1 and the base flow for the radially evolving liquid is simply obtained using the 
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incompressibility equation. Since the base flow depends only on the radial direction the 
pressure can be integrated using eq. (J1J, which gives: 

uoj = ~, v o = and p 0j = -r]jk 2 ln(r) + P 0j . 

The solution for the pressure field depends on the constant P i — P02 — 1/C°ff + 
which matches the jump in normal stress according to Laplace Law, due to surface tension 
at the interface. 

Assuming a perturbation expansion of the non-dimensional interface position p = po + 
ef(t)e' ine and = \& j + e ^ij> with a time dependent amplitude f(t) and normal mode 
Ansatz of wave number n. 

The Brinkman equation being linear, the general solution for eq. fl6]) is used. 

*n = i(o„^^ + t„r") e -" 9 , (10) 

The general solution is build of powers of r and modified Bessel functions of first and 
second kind, I n and K n . The parameters a n , b n , c n and d n have to be determined from 
the boundary conditions at the interface. The structure of eq. (EJ) shows that the classical 
solutions of the Laplace equation are also solutions of the Brinkman equation. Applying 
flattened boundary conditions at the interface that consider the first order Taylor expansions 
around the original interface position. 

Impermeability: The interface moves with the inner radial velocity such that p + 
ef e ind = uoi + eun, where the dot symbol denotes differentiation with respect to time. This 
yields 

fe M = ^fe M --^ + 0(e) , (12) 
or r 09 r=i 

which further leads to: 

n(a n + b n )-f-f = 0. (13) 

Continuity of radial velocity: Matching the outer velocity to the inner velocity at the 
interface, Un = w 12 : 

89 36 r =i 1 j 

leads to 

a n + K - c n - d n = 0. (15) 



Continuity of angular velocity: The equality of tangential velocities on the flattened 
interface, V\\ = vu- 

dr dr r=i 

leading to: 

, In+i(k) K n+ i(k) . 

a n k r n\ + c ™ k ts n\ + n \ an + bn ~ °n + d„) = 0. (17) 
l n (k) K n {k) 

Continuity of tangential stress: To ensure continuity of the tangential stress at the 
interface eq. (jSJ): 

/ d 2 ip d_(J_^ \\ 

\ r 2 d6 2 11 <9r \r<9r 

/ Q 2 d / d \\\ 

\ r 2 d6 2 12 ^ dr\rdr 12 //L=i ^ ^ 



Which simplifies to: 

77 ( 2a n n 2 + 26 n n 2 + a n k 2 - 2a n - 2a n n - 26 n 

\ I n [k) / 

— 2c n n 2 — 2c?„ n 2 — c n k 2 — 2c n k — n+1 ^ _|_ 2c n n — 2d n n — 0. (19) 

K n {k) 

Laplace Law: The jump of the normal stress balance is more involved since it requires 
to recover the pressure, which had been left out in the calculation. Integrating eq. flH]) for 9 
relates the pressure perturbations to the general solutions. 



Pa = w J o {- 



(n 2 , , 1 \ d^n 2n 2 \ . . 

which simplifies into 

Pll = - v b n r n k 2 e me , and p 12 = d n r- n k 2 e mB . 
The radial stress balance in order e for the flattened interface according to eq. (Q: 

or or \rov / or z 

ar ar \ray / or 1 
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Recalling that the factor 7r/4 in front of 1/Ca stems from a corrected boundary condition 
by Park and Homsy 17 . Inserting the expressions for velocity and pressure results in: 



+ b n k 2 — 2a n n — 2b n n + 2a n nk + 2a n n 2 

In{k) 



+2b n n 2 + 4/(i)) + d n k 2 - f{t)k 2 + 2c n n + 2d n n 

+ 2c n nfc i ^ + . 1 / ( ^ ) - 2c n n 2 + 2d n n 2 - 4f(t) = - n 2 ). (22) 

K n {K) 4Ga 

Using the boundary condition equations ( fl3l) . (115j) . (fl7j) . ( 1191) and (1221) the unknowns 
/, a n , b n , c n , d n have to solve an eigenvalue problem, with / as a Eigenvalue. Assuming that 
the growth of the instability is much faster than the evolution of the base flow, it is assured 
that the aspect ratio k does not change with time, which leads to the possibility to express 
the function f(t) = f exp(at). In this case f/f — cr, where a is the growth rate. 

IV. RESULTS GROWTH RATES FOR LINEAR THEORY 

One is interested in the wave number n with the highest growth rate a for given param- 
eters k and Ca. In figure [2] the most unstable wave numbers plotted in the same way as 
done by Maxworthy are shown, including the results from our stability analysis. Maxworthy 
chose to plot the results for a modified wave number being: 

~l2n r 




'max 



Ca k 

because then Paterson's result is approximately represented as a straight line. With the 
corrected boundary conditions of Park and Homsy this line should be constantly 2.26 for 
all capillary numbers. Maxworthy showed, as the capillary number approaches one the 
observed number of fingers becomes much smaller than predicted by Paterson. Our findings 
are represented by the full line and show an improved agreement. 

Maxworthy's data is spread, which can be explained by the relatively flat summit of the 
dispersion relation, as seen figure HI which does not cause a sharp wavelength selection. 
Displaying also possible fingers that have a growth rate of no less than 1% of the maximal 
growth rate. The spread takes an almost symmetric shape about the maximal growth rate. 

The curves vary slightly for low aspect ratio k < 10 2 . But in experimental applications 
the aspect ratio was chosen sufficiently large, k > 10 2 so the graph collapses to a single curve 
for all values of k in the range of interest. 




-2 -1.5 

log(Ca) 



FIG. 2. Most unstable modified wave number against logarithm of capillary number. Plotted 

results are: — • — Potential flow, Viscous potential flow and — Brinkman model with aspect 

ratio k = 200. Maxworthy's experimental results are also shown, indicated by o symbols. The grey 
shaded region includes all wave numbers of the Brinkman model, where the growth rate is within 
1 % of the maximal growth rate. 

The dispersion relation that is obtained from the Eigenvalue problem is for brevity given 
only for the case where the inner viscosity is zero (the general expressions are given in the 
appendix). The growth rate a becomes: 



A(k,n) + fCa- 1 J B(A;,n) 



(23) 



C(k,n) 

In this form it can be easily seen that the growth rate depends on a parameter pair, where 
A(k,n) weights the flow rate and B(k,n) weights the influence of surface tension. The 

10 




10" 



10 



10" 

Capillary number 



10 z 



FIG. 3. Most unstable finger number as a function of capillary number: Potential flow, 

— • — Viscous potential flow and — Brinkman model. Dotted lines show the experimental result of 
Maxworthy for low surface tension as deduced by Homsy 1 ^ and the theoretical result of Logvinov 
et al. who used the Brinkman model for zero surface tension fingering of a flat interface with aspect 
ratio k = 200. For higher capillary number the number of fingers becomes proportional to the gap 
width. 

parametric functions are: 

A(k, n) = k 4 n - k 4 - 8n 2 - 8n 4 + Ak 2 n - 8k 2 n 2 
K n+1 (k) 



+ 



K n (k) 



-2k 6 + 2k 6 n + 4fcn + Akrf 



B(k, n) = k 2 n - k 2 n 3 + An 4 - An 2 + Kn+ )[ k ) ' (2kn - 2kn^ 

K n (k) 



C(k, n) = k 4 - 8n 2 + 8n 4 + Ak 2 n 2 + Kn+ )\ k ) (2k 3 - Akn 3 + Akn). 

K n (k) 



(24) 
(25) 

(26) 



For a better comparison with existing results the growth rate a obtained by Paterson^ is 
decomposed in the same way. Correcting only the influence of the curvature in Paterson's 
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result by factor 7r/4, which is done for the corrected boundary condition. 

Reorganizing the equation in the form of eq. ( l23i) gives: 

A(k,n) = (n - l)k\ B(k,n) = n(l - n 2 )k 2 
and C{k, n) = fc 4 . 

80- 
70 

60 - 




Wave number n 

FIG. 4. Dispersion relations using different models: • • • Potential flow, Viscous potential flow 

and — Brinkman model. The aspect ratio is fixed to k = 400, while the capillary number changes 
between 0.032,0.1,0.32 and 1.0. 

This rephrasing shows that the results of Wilson 4 - and Paterson^ are recovered in the limit 
of large aspect ratio k. However it is apparent that for zero surface tension, where Paterson's 
pairing of A(k, n)/C{k, n) grows linearly with n, our relation has stabilizing terms that grow 
faster in wave number n than the destabilizing terms. This ensures that for finite aspect 
ratio k, there will always be a finite cut-off wave number n. 
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Plotting in figure H] the growth rate a as a function of the wave number for different 
capillary numbers at constant aspect ratio k = 400 shows two properties in the dispersion 
relation induced by the Brinkman model: The growth rates are smaller for all wavelengths 
and the maximum growth rates shift to lower wave numbers as higher modes receive larger 
damping. 

The number of fingers for zero surface tension has shown to depend only on the gap 
width or aspect ratio. Thus plotting the number of most unstable finger divided by the 
aspect ratio k against capillary number in figure [3j For k > 10 2 the results approximately 
collapse on a single curve for all values of k. For large capillary numbers the prediction tends 
asymptotically to the result of Logvinov et al.- 2 who studied zero surface tension flows in 
rectangular channels using the Brinkman equation. The experimental result, identified from 
Maxworthy's data, is missed by a factor 2. Not surprisingly the result of Paterson gives a 
far too high number of fingers. 

V. CONCLUSION 

A dispersion relation for Saffman- Taylor modeled by the Brinkman equation covers a 
wide range of capillary numbers and shows a good approximation with experiments. The 
decreased wave number of the maximal growth is achieved by keeping the in-plane stresses. 
When the surface tension has only little influence, the in-plane stresses decrease the growth 
rate for high wave numbers. The high wave numbers create large velocity gradients and 
thus shear stresses at the interface. Darcy's law can not account for them, in contrast to the 
Brinkman model, which allows for a cut-off wave number even with zero surface tension. 

The finger size for zero surface tension is found to scale like the cell thickness h, as 
found by Logvinov et al.— who considered a rectangular cell and by Kim et ai— for radial 
injection. In the contrast to the latter study, based on viscous potential flow, our solution 
satisfies continuity of tangential velocity and shear stress perturbations. 

The mismatch that is still observed at low capillary numbers might be attributed to 
the formation of thin films on the advancing frontr^. On one side, the deposition of films 
reduces the gap width and speeds up the fluid. On the other side, the Laplace law should 
be modified to account for the strong curvature change in the dynamical meniscus region. 
Whether including these effects does improve the predicted finger number selection is the 
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subject of our current work. 

Another natural perspective is to consider the full nonlinear solution of the radial fingering 
by the numerical resolution of the two-phase Brinkman equations. While classical tip- 
splitting phenomena are recovered, our preliminary results point to a strong dependence of 
the wavelength selection on the initial conditions and background noise, a large statistical 
deviation is being observed. 
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APPENDIX: FULL DISPERSION RELATION 

The growth rate for wave number n depending k, Ca and n. 

a = -(-(4/x k 4 Ca K ) - (8h k 3 Ca K x ) - (32h Ca n 2 K ) - 4I^n 2 K Q + 4I l7 rn 4 K 
-(32/i Can 4 K ) + (I6/1 k 2 CanK ) - (4hk 4 r]CaK ) - (32 1 1 k 2 Can 2 K ) 
+ (4/i k 4 n Ca K ) + (8/1 k 3 nCaK 1 ) + (16h knCaK l ) + I x k 2 nirK + 2I X k wkK x 
-hk\n 3 K -2hkTTn 3 K x + (I6/1 kCa n 3 - (4k 4 I nCa + (8/1 k 3 rj 2 Ca Ki) 
-(4k 4 I r] 2 Ca Kx) + (64/^n 2 K Ca) + (64/^n 4 Ca K ) - (32I l i 1 2 n 2 K Ca) 
+4rihirn 2 K - 4i 1 I l im 4 K - (32n 2 n 4 Ca h K ) - (32/x k nnCa K x ) + (48/! k 2 r]Ca n 2 K ) 
-(4I X k 4 nnCa K ) - (16/! A; 3 nnCa K x ) - (32^ knn 3 Ca - (32n 3 I Q nk Ca K ) 
-(8n 2 I r]k 3 Ca K ) + (16n 2 I nk 2 Ca K x ) - (8n I nk 3 Ca K ) - (lQn 2 I n 2 k 2 Ca K x ) 
+ (32n 3 1 r] 2 kCaK ) + (4k 4 I G nnCaK l ) + (\U x knn 2 CaK x ) - (I6/1 k 2 rfCan 2 K ) 
+ (8hk 3 nr] 2 Ca K x ) + (8n 2 I n 2 k 3 Ca K ) + (8n I G n 2 k 3 Ca K ) - (4k 4 I m?Ca K ± ) 
-(I6/1 k 2 rfCa n K ) - 2rjh k mxK x + 2nhknn 3 K x - 2nn 2 I nk K + 2nn 4 1 nk K 

+nk 2 I nir K x - nk 2 I im 3 K x + (16^ 2 n 3 Ca h k K 1 ))/(4Ca) 
/(-hk 4 r]K - 4hk 2 n 2 K - 4I x knK x + 4I 1 kn 3 K 1 - k 4 hnK x + 2I l k 3 r ] 2 K 1 - k 4 Itf] 2 K x 
-16hrin 2 K + lQhr]n 4 K + 8hn 2 n 2 K - 8rfn 4 hK - hk 4 K - 2I 1 k 3 K 1 + 8hn 2 K 
-8hn 4 K - 4nI r]K l k 2 + 8I l knnK l + 4I 1 k 2 nn 2 K - 8I 1 knn 3 K l - 8n 3 I r]kK + 4n 2 I nk 2 K l 
-2nI r]k 3 Ko - 4n 2 I n 2 k 2 K 1 + 8n 3 I ri 2 kK - 4I 1 knrj 2 K l + 2nI ri 2 k 3 K - 4I x k 2 n 2 nK^ 
+4rj 2 n 3 I l kK 1 + 8n 2 I r]kK - 8n 2 I r] 2 kK + 4nI r ] 2 K 1 k 2 + 4hk 2 nnK ). 
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